%% 逆牛顿类算法Quasi Newtonian algorithm
%% 参数初始化
lb=[-1e2,-1e2];ub=[1e2,1e2];%变量范围
n=length(lb);% 变量个数
alpha0m=1/n*mean(ub-lb);% 初始步长
rand1=20*rand()+0.1;
f=@(x)(x(1)-1)^2+(x(2)-2)^2;% 目标函数
x_min=[1,2];% 最优解
x0=rand1*x_min;% 令迭代初始点在最优解附近
%% 牛顿类算法搜索计算过程
alpha0=alpha0m;% 初始步长
mn=300;% 迭代次数
e=1e-3;% 判断梯度为0点的阈值
gama2=1/mn;% 迭代次数倒数
mn1=100;% 非单调准则判断迭代次数
x=x0;% 初始点位置
dx=(1e-6)*(ub-lb);% 求梯度时维度微元设置
gama1=0.9;% 回退法参数gama1设置
c1=1e-3;% 常数c1设置
c2=1e-3;% 常数c2设置
alpha=alpha0;% 迭代步长初始化
for i=1:mn
    alpha0=(mn-i+1)*gama2*alpha0m;% 该步迭代最大步长限制
    alpha=alpha0;
    tic
    c_u=zeros(1,n);% 该点梯度初始化
    for bianliang=1:n% 求梯度
        c_un1=x(bianliang)-dx(bianliang);
        c_un2=x(bianliang)+dx(bianliang);
        if x(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
            % 量下限时处理
            r_un=rand();
            c_un1=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
        elseif x(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
            % 于变量下限时处理
            r_un=rand();
            c_un2=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
        end
        x1=x;x2=x;
        x1(bianliang)=c_un1;x2(bianliang)=c_un2;
        c_u(bianliang)=(f(x2)-f(x1))/(2*dx(bianliang)+eps);
    end
    G1=c_u;% 该点梯度
    if(norm(c_u)<e)% 当该点处梯度足够小时，说明找到函数驻点，跳出循环
        break
    end
    c_u=c_u/(norm(c_u)+eps);% 确定该点最终梯度方向，当计算中未出现超过变量限制时为
    % 梯度方向，否则为次梯度方向。
    for j=1:mn1
        if i==1
            dk=-c_u;
            xk=x+(dk*alpha);% 迭代点
        else
            dk=(-(B+eps)\c_u')';
            dk=dk/(norm(dk)+eps);
            if(isnan(dk(1)))
                dk=zeros(1,n);
            end
            xk=x+(dk*alpha);
        end
        for bianliang=1:n% 该步限制迭代之后的点位于自变量范围之内
            if xk(bianliang)<lb(bianliang)
                r_un=rand();
                xk(bianliang)=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
            elseif xk(bianliang)>ub(bianliang)
                r_un=rand();
                xk(bianliang)=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
            end
        end
        c_u1=zeros(1,n);% 该迭代点梯度初始化
        for bianliang=1:n% 求该迭代点梯度
            c_un1=xk(bianliang)-dx(bianliang);
            c_un2=xk(bianliang)+dx(bianliang);
            if xk(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
                % 量下限时处理
                r_un=rand();
                c_un1=r_un*lb(bianliang)+(1-r_un)*xk(bianliang);
            elseif xk(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
                % 于变量下限时处理
                r_un=rand();
                c_un2=r_un*ub(bianliang)+(1-r_un)*xk(bianliang);
            end
            x1=xk;x2=xk;
            x1(bianliang)=c_un1;x2(bianliang)=c_un2;
            c_u1(bianliang)=(f(x2)-f(x1))/(2*dx(bianliang)+eps);
        end
        G2=c_u1;% 梯度
        % 该步添加Wolfe搜索准则
        if (f(xk)<=f(x)+c1*alpha*(c_u*dk')&&G2*dk'>=c2*G1*dk')% Wolfe准则
            break
        else
            alpha=gama1*alpha;
        end
    end
    c_u1=zeros(1,n);% 迭代点梯度初始化
    for bianliang=1:n% 求迭代点梯度
        c_un1=xk(bianliang)-dx(bianliang);
        c_un2=xk(bianliang)+dx(bianliang);
        if xk(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
            % 量下限时处理
            r_un=rand();
            c_un1=r_un*lb(bianliang)+(1-r_un)*xk(bianliang);
        elseif xk(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
            % 于变量下限时处理
            r_un=rand();
            c_un2=r_un*ub(bianliang)+(1-r_un)*xk(bianliang);
        end
        x1=xk;x2=xk;
        x1(bianliang)=c_un1;x2(bianliang)=c_un2;
        c_u1(bianliang)=(f(x2)-f(x1))/(2*dx(bianliang)+eps);
    end
    G2=c_u1;% 梯度
    % 利用BB算法计算步长
    sk=xk-x;yk=G2-G1;
    if(i==1||det(B)==0)
        B=(yk'*yk)/(sk*yk'+eps);
    else
        B=B+(yk'*yk)/(sk*yk'+eps)-((B*sk')*(B*sk')')/(sk*B*sk'+eps);
    end
    if(isnan(xk(1)))
        break
    end
    x=xk;
    toc
end
f_min=f(x);